#DESCRIPTION:
#Kinematic analysis of geologic structures vis-a-vis slopes


#1. libraries and modules

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import beta,gamma,gaussian_kde, kstest, ks_2samp, ttest_ind,norm
import matplotlib.cm as cm

import BartonBandis as BB

##
##def generate_sample_grad_asp_data(num_data):
##    grad=beta.rvs(2,5,loc=0,scale=75,size=num_data)
##    asp1=beta.rvs(2,5,loc=-20,scale=220,size=num_data/2)
##    asp2=beta.rvs(2,5,loc=160,scale=220,size=num_data/2)
##    asp=np.concatenate((asp1,asp2))
##    asp=np.where(asp<0,360+asp,asp)
##    asp=np.where(asp>=360,asp-360,asp)
##    return np.round(grad,1),np.round(asp,1)

def app_dip(df,d_dip, d_ddir):
    delta_az=np.abs(df.asp.values-d_ddir)
    delta_az=np.where(delta_az>180,360-delta_az,delta_az)
    d_adip=np.where(delta_az<=90,
                    np.degrees(np.arctan(np.tan(np.radians(d_dip))*np.cos(np.radians(delta_az)))),
                    -1*np.ones(delta_az.shape))
    return np.round(delta_az,1), np.round(d_adip,1)

def KinematicAnalysis(df,d_dip, d_ddir):
    delta_az, d_adip=app_dip(df,d_dip, d_ddir)
    KA=np.where((d_dip>=0)*(d_dip<df.grad.values)*(delta_az<=20),
                np.ones(d_adip.shape),
                np.zeros(d_adip.shape))
    df['delta_az']=delta_az
    df['d_adip']=d_adip
    df['KA']=KA
    return df

##def spherical_to_polar_to_cart(dip,ddir):
##    #spherical
##    dip=np.asarray(dip)
##    ddir=np.asarray(ddir)
##    z=1+np.sin(np.radians(dip))
##    x=np.cos(np.radians(dip))
##    zenith_angle=(np.arctan(x/z))
##    #to polar
##    r=np.tan(zenith_angle)#/(1-np.cos(zenith_angle))
##    theta=-ddir+90
##    #to cartesian
##    X=r*np.cos(np.radians(theta))
##    Y=r*np.sin(np.radians(theta))
##    return X,Y
##
##def plot_stereonet():
##    ax=plt.gca()
##    theta=np.linspace(0,2*np.radians(360),1000)
##    X=np.cos(theta)
##    Y=np.sin(theta)
##    ax.plot(X,Y,'k-')
##    ax.plot([-1,1],[0,0],'k-',lw=0.1)
##    ax.plot([0,0],[-1,1],'k-',lw=0.1)
##    for d in [15,30,45,60,75]:
##        z=1+np.sin(np.radians(d))
##        x=np.cos(np.radians(d))
##        zenith_angle=(np.arctan(x/z))
##        r=np.tan(zenith_angle)
##        X=r*np.cos(theta)
##        Y=r*np.sin(theta)
##        plt.plot(X,Y,'k-',lw=0.1)
##    ax.axis('equal')
##    ax.set_xticks([])
##    ax.set_yticks([])
##
##def draw_great_circle(dip,ddir):
##    a=np.linspace(-90,90,1000)
##    adipdir_list=ddir+a
##
##    delta_az=np.abs(adipdir_list-ddir)
##    delta_az=np.where(delta_az>180,360-delta_az,delta_az)
##    d_adip=np.where(delta_az<=90,
##                    np.degrees(np.arctan(np.tan(np.radians(d_dip))*np.cos(np.radians(delta_az)))),
##                    -1*np.ones(delta_az.shape))
##
##    X,Y=spherical_to_polar_to_cart(d_adip,adipdir_list)
##    ax=plt.gca()
##    ax.plot(X,Y,'r-')
##    
##    


##
##
###generate sample gradient data
##num_grad_data=5000
##grad,asp=generate_sample_grad_asp_data(num_grad_data)
##
##path_to_file='/home/eglsais/Dropbox/Thesis/csv_files_originals/'
##df_in=pd.read_csv(path_to_file+'DF_el_li_gr_asp')
##
##
##grad=df_in.grad.values
##asp=df_in.asp.values
##
##print grad.shape, asp.shape
##
##
##d_dip_list=[40,50]
##num_disc=len(d_dip_list)
##fig,ax=plt.subplots(nrows=2,ncols=num_disc, sharex=True,sharey=True)
##for d in range(len(d_dip_list)):
##
##    #sample discontinuity data
##    d_dip=d_dip_list[d]
##    d_ddir=60
##
##    #create dataframe
##    df=pd.DataFrame()
##    df['s_grad']=grad
##    df['s_asp']=asp
##
##    #conduct kinematic analysis for given discontinuity
##    df=KinematicAnalysis(df,d_dip, d_ddir)
##
##    #selecting rows for FS analysis
##    df_FS=df[(df.KA==1)]
##    print df_FS
##
##    ####################################
##    #conduct BB FS analysis
##
##    #defining constants
##    uw=24*10**-3 #MN/m3
##    base=30
##    Ln=5
##
##    #defining slope gradient (beta) and discontinuity dip (alpha)
##    beta=df_FS.s_grad.values
##    alpha=df_FS.d_adip.values
##
##    #defining random variables
##
##    numrand=1000
##    sign_min,sign_max=0.,3.     #MPa
##    phir_min,phir_max=10,25     #degrees
##    JRC0_min,JRC0_max=0,20     
##    JCS0_min,JCS0_max=10,60     #MPa
##
##    sig_n,phi_r,JRC0,JCS0=BB.generate_random_BB_vars(sign_min,sign_max,
##                                                  phir_min,phir_max,
##                                                  JRC0_min, JRC0_max,
##                                                  JCS0_min,JCS0_max,
##                                                  numrand)
##
##    df_BBFS=BB.BB_stability_analysis(uw,base,Ln,
##                              beta,alpha,
##                              phi_r,JRC0,JCS0,0)
##
##    df_FS['FS_mean']=df_BBFS.FS_mean.values
##
##    df_FS['FS_std']=df_BBFS.FS_std.values
##
##
##    df_FS['FS_pcF']=df_BBFS.FS_pcF.values
##
##    print df_FS
##    ############################################################
##
##    #plotting KA results
##
##
##
##    
##    #plotting all slope data
##    X1,Y1=spherical_to_polar_to_cart(df.s_grad.values,df.s_asp.values)
##    plt.sca(ax[0,d])
##    plot_stereonet()
##    ax[0,d].scatter(X1,Y1,s=1,marker='.')
##    draw_great_circle(d_dip,d_ddir)
##
##    #plotting kinematically admissible slope data
##    X2,Y2=spherical_to_polar_to_cart(df_FS.s_grad.values,df_FS.s_asp.values)
##    
##
##    plt.sca(ax[1,d])
##    plot_stereonet()
##    #ax[1,d].scatter(X1,Y1,s=1,c='0.6',marker='.')
##    draw_great_circle(d_dip,d_ddir)
##
##    col=df_FS.FS_pcF.values
##    g=ax[1,d].scatter(X2[np.argsort(col)],Y2[np.argsort(col)],c=np.sort(col),cmap=cm.jet, vmin=0,vmax=100,linewidth=0)
##    #plt.colorbar(g)
##
##fig.tight_layout()
##
##plt.show()
##
##
##



        
    



    
